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We present an algorithmic study for the simulation of two massless flavors of 0(a) improved Wilson quarks with 
Schrodinger functional boundary conditions. The algorithm used is Hybrid Monte Carlo with two pseudo-fermion 
fields as proposed by M. Hasenbusch. A gain in CPU cost of a factor two is reached when compared to one 
pseudo-fermion field due to the larger possible step-size. This study is integrated in the ALPHA project for the 
computation of the running of the renormalized quark mass. We include an update on these physics results. 



1. HMC with two pseudo-fermion fields 

We study a variant of the Hybrid Monte 
Carlo (HMC) algorithm that uses two pseudo- 
fermion fields per degenerate flavor doublet, as 
proposed by M. Hasenbusch 2^3; and recently 
tested in |l . The Wilson-Dirac operator with 
Schrodinger functional (SF) boundary conditions 
Dyj + rriQ is considered together with 0(a) im- 
provement 5D respectively T (clover term) and 
even- odd preconditioning 
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For one mass-degenerate flavor doublet of quarks 
the partition function reads 



Z = 



exp(-S'g(C/)) dct(l + Tee)^ dot Q2 ^ (3) 



Q = 5075{1 + Too - Moe{l + TeeV^Meo} , (4) 

with CQ = {1 + 64k2)-i and 5g(;7) is the Wilson 
plaquette action. To simulate cq. JSJ we use the 
Hybrid Monte Carlo algorithm. The determinant 
det is represented in terms of pseudo-fermion 
flelds (p. If Q is factorized in npf factors then 
one pseudo-fermion fleld for each factor is intro- 
duced. The aim of the factorization is to make 
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the fermionic contribution to the forces associ- 
ated with each factor as small as possible. We 
take npf — 2 and the factorization 

Q = {Q){Q-^Q). Q = Q-ip, (5) 

which leads to the pseudo-fermion actions i — 1,2 
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The real parameter p is chosen to minimize the 
sum of the condition numbers of the two opera- 
tors appearing in eq. ©, 



(Q')A,„ax(Q') 



1/4 



(8) 



in terms of the smallest and largest eigenvalue 
of Q^. The operators in eq. © have then both 
condition number equal to the square root of the 
condition number of Q^. 

This algorithmic study is part of the ALPHA 
large scale simulations for the computation of the 
running of the renormalized quark mass For 
the temporal and spatial lattice size, the back- 
ground gauge fleld and the parameter controlling 
the spatial boundary conditions of the fermion 
fields we take respectively 



T = L, C = C' = 0, 



= 0.5. 
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Figure 1. Step-sizes St for 80% acceptance. Open 
symbols are for Upf = 1, filled symbols and crosses 
for Tipf = 2. 

The theory is simulated along the critical line 
where the PC AC mass vanishes. These simula- 
tions are possible in the SF due to the infrared 
cut-off proportional to in the spectrum of 

the Dirac operator squared. In our simulations 
the renormalized coupling u ~ g^{L) [H] takes 
values in the range u — 1.0 . . . 5.7. In the SF u is 
presumably a monotonically growing function of 
L and our simulations correspond approximately 
to the range L — 10"^ ... 1 fm. 

2. PERFORMANCE 

In the integration of the molecular dynamics 
equations of motion we adjust the step-size St 
to yield an acceptance of 80% for trajectories of 
length one. Fig.^shows the inverse step-size, i.e. 
the number of steps, as a function of the average 
condition number of 



{k) = 
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Data are shown for different lattice sizes L/a 
and renormalized couplings u = 1.0,1.1 (trian- 
gles), u = 2.5 (squares), u = 3.3 (pentagons), 
u — 5.7 (circles) together with simulations at 
(3 ~ 5.2, K = 0.1355 (crosses). The filled sym- 
bols and crosses are obtained with ?ipf = 2 and 
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Figure 2. Integrated autocorrelation time of Zp. 
Data are for lattices L/a = 12 (triangles), L/a = 
16 (squares) and L/a = 24 (circles), open symbols 
for ripf = 1 and filled symbols for Upf = 2. 

the open symbols with Upf = 1 (standard HMC). 
With Upf = 2 we can choose step-sizes which are 
a factor two larger than with ripf = 1. 

In Fig. 13 we plot the integrated autocorrelation 
time Tint(^p) of the pseudoscalar density Zp 
as a function of u. Zp is the quantity we need 
in order to compute the running of the renor- 
malized mass. The units are CPU seconds per 
lattice point for the simulation of one lattice on 
one APEmille crate, which consists of 128 nodes 
and has a peak performance of 68 GFlops. Data 
are for lattices L/a = 12 (triangles), L/a = 16 
(squares) and L/a = 24 (circles), open symbols 
for Upf = 1 and filled symbols for ripf = 2. A sys- 
tematic difference between ripf = 1 and Upf = 2 is 
seen for the L/a — 16 lattices, the smallest cou- 
plings show a reduction in CPU cost by about a 
factor two for Upf = 2. This CPU gain is due 
to the larger possible step-size. A comparison at 
L/a — 24 would demand too much CPU time, 
precisely because of this it was essential to our 
project to speed up the standard HMC. 

3. THE RUNNING MASS 

The running of the renormalized quark mass 
rn{p,) with the renormalization scale /i = 1/i is 
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Figure 3. The running of the renorniahzed quark 
mass. 

extracted from the step scahng function cr(u) of 
the couphng and the step scahng function crp(u) 
of the pseudoscalar density, 
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CTp(w) 
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We solve the joint recursion: fc = 0, 1, 2, . 
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which gives the running of rn{^) starting at the 
scale 1/Lmax defined by 5^(Lmax) = 4.6800. 
For the step scaling function crp (u) we take the 

average Sp^' (u) of its values Sp(w, a/L), which 
we computed on L/a = 8, 12 (and 2L/a) lattices 
for 6 couplings in the range m = 0.98...3.3. Then 
we interpolate by the Ansatz 



Ep'iu) = l-lii{2)dou+p2u'^ +P3U^ (18) 



with do perturbative, fitted parameters P2 and pa 
and use this fit formula in the recursion eq. (|16|l . 
At high energies contact with the perturbative 
regime is made from which the RGI parameters 
A and M can be extracted [l] 

In(Ai^ax) = -1.34(7) (19) 

isra = '■="'<='>■ <^»' 

Fig. 13 is then obtained from the coefficients 
eq. ((T3)l and eq. l(T7|) . The errors of the points 
come from the statistical errors in the coefficients 
Wk , the scale ambiguities in the quantities ALmax 
and M/TO(1/Lniax) are not shown. 

The outlook of this work is the determination 
of imax uslug available data from JLQCD |H] and 
UKQCD A hadronic scheme has to be set up 
to determine combinations of M„, Md and Mg. 
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